R Notebook

Libraries

Load data

I am using the AHA dataset for writing my code because it is the simplest dataset I can work with. Simplest = significant number of baseline samples without issue of repeated measures, balanced distribution of classes, likely difference to be detected present.

p_ps = readRDS("/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/POMMS/R24_POMMS/data/processed/20221104_pomms_metadata_plant.rds")

p_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 252 taxa and 384 samples ]
sample_data() Sample Data:       [ 384 samples by 71 sample variables ]
tax_table()   Taxonomy Table:    [ 252 taxa by 8 taxonomic ranks ]

sample_data(p_ps)

Sample Data:        [384 samples by 71 sample variables]:

Preprocess

Remove problematic samples

Taking out samples with less than or equal to 0 reads and keeping only entry timepoint reduces samples from 384 down to 240.

p_ps = p_ps %>%
  subset_samples( reads > 0) %>%
  subset_samples(timepoint == "Entry")

p_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 252 taxa and 240 samples ]
sample_data() Sample Data:       [ 240 samples by 71 sample variables ]
tax_table()   Taxonomy Table:    [ 252 taxa by 8 taxonomic ranks ]

Class distribution in this dataset: Case = 193 Control = 47

sample_data(p_ps) %>%
  as.data.frame() %>%
  as_tibble() %>%
  group_by(group) %>%
  summarise(count = n())

Data transform

otu_clr = abundances(p_ps, "clr") %>% # rclr transform could not be performed, get Nans
  t()

p_ps = phyloseq(otu_table(otu_clr, taxa_are_rows = FALSE),
                   sample_data(sample_data(p_ps)),
                   tax_table(tax_table(p_ps)))

rm(otu_clr)

Feature table

Extract data

Order of asv colnames in OTU table and tax table is the same.

features = p_ps@otu_table %>%
  as.data.frame() %>%
  as_tibble(rownames = NA) %>%
  rownames_to_column("sample")

feature_labels = c("sample")

 names = tax_table(p_ps) %>%
      data.frame() %>%
      lowest_level() %>%
   pull(name)

 feature_labels = append(feature_labels, names)

 feature_labels = make.unique(feature_labels) %>%
   make.names()

 colnames(features) = feature_labels

features

Attach data labels

sample_labels = p_ps@sam_data%>%
  as.data.frame() %>%
  as_tibble(rownames = NA) %>%
  rownames_to_column("sample") %>%
  select(sample, group)

features = sample_labels %>%
  left_join(features)

Joining, by = "sample"

Zero variance features

No features have zero variance

nzv = features %>%
  select(!c(sample,group)) %>%
  nearZeroVar(saveMetrics = TRUE)

rm(nzv)

Correlated predictors

feature_cor = features %>%
  select(!c(sample,group)) %>%
 cor() 

summary(feature_cor[upper.tri(feature_cor)])

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.43611 -0.03523  0.06750  0.13015  0.22513  1.00000 

high_cor = findCorrelation(feature_cor, cutoff = .75)

filtered_features = features[,-high_cor]

rm(high_cor)
rm(feature_cor)

Linear dependencies

feature_ld = filtered_features %>%
  select(!c(sample,group)) %>%
  findLinearCombos()

filtered_features = filtered_features[, - feature_ld$remove]

Input table

Preprocessing reduces column number to 131 from 235

input = filtered_features %>%
  mutate(group_shuffle = sample(group) %>%
           as.factor(),
         group = as.factor(group))

Modeling

Load functions from classification_function_AA file

results = loop_label(iterations =20, input, label_list, percent = 0.7, downsampling = TRUE)

[1] "will downsample on training data"

Warning: There were missing values in resampled performance measures.

[1] "group"
[[1]]
[1] "group"

[[2]]
          Reference
Prediction Case Healthy
   Case     676     125
   Healthy  464     155

[1] "will downsample on training data"

Warning: There were missing values in resampled performance measures.Warning: There were missing values in resampled performance measures.

[1] "group_shuffle"
[[1]]
[1] "group_shuffle"

[[2]]
          Reference
Prediction Case Healthy
   Case     578     142
   Healthy  562     138

rm(label_list)

Plotting

AUC plot and t-test

auc_df %>%
  pivot_longer(cols = c("group", "group_shuffle"), names_to = "feature", values_to = "AUC") %>%
  ggplot()+
  geom_boxplot(aes(x=feature, y = AUC))+
  theme_classic()

t.test(auc_df$group, auc_df$group_shuffle)

    Welch Two Sample t-test

data:  auc_df$group and auc_df$group_shuffle
t = 3.6882, df = 37.743, p-value = 0.0007082
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.04036328 0.13863672
sample estimates:
mean of x mean of y 
   0.5815    0.4920 

Importance plots

look-up table

#Making asv reference for names in model
lookup_asv = tax_table(p_ps) %>% 
     data.frame() %>%
      lowest_level() %>%
   rownames_to_column("asv") %>% 
   select(asv, name) %>%
   mutate(name = make.names(make.unique(name)))

features

summarize importances

top 10 ranked list

top10

 [1] "Helianthus.annuus" "Brassica"          "Solanaceae"       
 [4] "Apiaceae"          "Capsicum.annuum"   "Poaceae"          
 [7] "NA."               "NA.3"              "Lactuca.sativa"   
[10] "Zea.mays"         

Importance plot

importance_df %>%
  pivot_longer(cols = c("NA.":"Apium.graveolens"), names_to = "food", values_to = "importance")  %>%
  filter(variable == "group" & food %in% top10) %>%
  ggplot()+
  geom_boxplot(aes(x= reorder(food,-importance), y = importance))+
  labs(x= "Food", y ="Importance") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust=1))

Some of the important NAs for predictions only map to bacterial species in BLAST.

Important taxa directions

lookup_asv %>%
  filter(name == "NA.") %>%
  pull(asv)

[1] "AAATCCTGTTTTATGAAAATAAACAAGGGTTTCATAAACCGAAAATAAAAAAG"